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Abstract: In this paper, the problem of estimating the attitude of a micro-nano satellite, 
obtaining geomagnetic field measurements via a three-axis magnetometer and obtaining 
angle rate via gyro, is considered. For this application, a QMRPF-UKF master-slave 
filtering method is proposed, which uses the QMRPF and UKF algorithms to estimate the 
rotation quaternion and the gyro bias parameters, respectively. The computational 
complexicity related to the particle filtering technique is eliminated by introducing a 
multiresolution approach that permits a significant reduction in the number of particles. 
This renders QMRPF-UKF master- slave filter computationally efficient and enables its 
implementation with a remarkably small number of particles. Simulation results by using 
QMRPF-UKF are given, which demonstrate the validity of the QMRPF-UKF nonlinear 
filter. 

Keywords: attitude determination; gyro; magnetometer; QMRPF-UKF; master-slave 
filtering 
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1. Introduction 

Spacecraft attitude is important for attitude control in many space missions. In the last four decades, 
a great number of research works have been devoted to the problem of estimating spacecraft attitude 
based on a sequence of noisy vector observations, resolved in the body-fixed coordinate system and in 
a reference system [1-7]. 

For the attitude sensor, the sun sensor, Earth sensor, star tracker, three-axis magnetometer, and 
gyro etc, are often used. The attitude estimation problem possesses an undesirable strong nonlinearity. 
The filtering-based methods that were developed in the 1980s embedded the attitude determination 
problem in the framework of stochastic filtering. The mostly used filtering method is the extended 
Kalman filter (EKF) [2,8]. Nevertheless, poor performance or even divergence arising from the 
linearization implicit in the EKF has led to the development of other filters [3,9,10]. Another often 
used method is the unscented Kalman filter (UKF), also known as the sigma-point filter [4]. Unscented 
filters are essentially based on second or higher-order approximations of nonlinear functions, which are 
used to estimate the mean and covariance of the state vector. In [4], an unconstrained three-component 
vector is used, based on the generalized Rodrigues parameters to represent an attitude error quaternion. 
The state vector includes the attitude error and bias vectors. UKF, as a Kalman filter mechanization, is 
sensitive to the statistical distribution of the stochastic processes driving the dynamic model: 
non-Gaussian distributions guarantee nonoptimality of the estimates. Recently, a new method using the 
particle filtering (PF) technique has been proposed for a spacecraft's attitude estimation [1,5-7]. Based 
on the concept of sequential importance sampling and the use of Bayesian theory, particle filtering is 
particularly useful in dealing with nonlinear and non-Gaussian problems. However, the most notorious 
disadvantage of a particle filter is its formidable computational complexity, since hundreds even 
thousands of particles are usually needed to achieve required approximation accuracy, and then it is 
difficult to achieve real time performance. A particle filtering algorithm [1] copes with the curse of 
computational complexity related to PFs by using an efficient initialization procedure, along with an 
importance weight cooling schedule and particle set portioning. This method is applicable to gyroless 
attitude determination setting, and the state vector is composed of a unit norm quaternion and a 
three-element angular-rate vector. In [5], particle filtering for attitude estimation using a minimal 
local-error representation (MLERPF) is proposed. The main character of MLERPF is that it does not 
work well when the number of particles is not large enough, which brings limitation to some 
applications. 

This paper develops a QMRPF-UKF master-slave filter and gives its implementation in the case of a 
three-axis stabilized micro-nano satellite, obtaining noisy geomagnetic field measurements via a 
three-axis magnetometer and obtaining the noisy angle rate via gyro. In order to reduce the 
computational cost rising from the increased state dimensions, the quaternion multiresolution particle 
filter (QMRPF) algorithm is used as the master filter to estimate the attitude quaternion, and the 
unscented Kalman filter algorithm is used as the slave filter to determine the gyro bias parameters. On 
the other hand, the new estimator copes with the curse of computational complexity related to the 
particle filtering technique by introducing multiresolution approach that permits a significant reduction 
in the number of particles [11]. This renders the QMRPF-UKF master- slave filter computationally 
efficient and significant computational savings of QMRPF master filter is achieved by drastically 
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reducing the needed number of particles. The performance of QMRPF-UKF in various scenarios is 
studied and compared with MLERPF method, which demonstrates the validity of QMRPF-UKF 
nonlinear filter. 

This paper is organized as follows. In Section 2, the system model of the attitude determination of a 
three-axis stabilized micro-nano satellite by using gyro and magnetometer is given. In Section 3, the 
QMRPF-UKF master-slave filter is described in detail. In Section 4, the QMRPF-UKF master-slave 
filter for the attitude determination problem is verified on the basis of its estimation accuracy and 
computational time. We draw some conclusions in Section 5. 

2. System Model 

2.1. Quaternion Process Model 

It is known that the body angular motion can be described in terms of the attitude quaternion by the 
following equation: 

q k+1 =n(G) k )q k (1) 

where <a k = \jo xk co yk <$ zk \ denotes the angular velocity vector of the rotation of 6 with respect 
to R . q k denotes the quaternion of rotation from a given reference frame R onto the body frame B at 
times k = 1, 2, . . . , oo . The quaternion is constructed from vector part q k and scalar part q Ak : 

q k =\.q k T q ik T (2) 

Assuming that co k is constant during the sampling time interval A?, the orthogonal transition matrix 
fl(a> k ) is expressed using co k and computed as follows: 

cos(0.5||«jAf)/ 3x3 -[^x] y/ k 

~¥k cos(0.5||<x>JA?) 



Q(co k ) = 



(3) 

4x4 



sin (0.511^ II At)co k 
Wk = — \T1\ — (4) 

IK II 

where [^x] denotes the cross-product matrix associated with the vector^ . 

2.2. Rate Sensor Measurement Model 

A common sensor that measures the angular rate is a rate-integrating gyro. For this sensor, a widely 
used model is given by: 

G>k=°>k + Pk +T lv,k (5) 

A + i = A + i u ,k (6) 

where G) k is the measured rate; /3 k is the gyro bias vector; T] v k and rj u k are independent zero-mean 
Gaussian white-noise processes with: 
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E{r lv (t)r J l(T)}=a^(t-r)I 3x3 (7) 
E{r lu {t)r 1 T u {T)}=a 2 u 8{t-T)I^ (g) 
where £{ • } denotes expectation and 8(t - r) is the Dirac_delta function. 

2.3. Observation Model 

A widely used attitude measurement model is given by: 

*>k= A (q k X+Sb k (9) 

where b k is the body-frame vector and r k is the reference-frame vector; 8b k denotes the measurement 
noise process, with known probability density function (PDF) denoted as Sb k ~ p Sbt ; where A(q k ) 

denotes the rotation matrix that brings the axes of R onto the axes of B at time k . The attitude matrix 
is related to the quaternion by: 

A(q k ) = [(q 4k ) 2 -q k T q k ]I 3x3 + 2q k q k T -2q 4k [q k x] (10) 



3. QMRPF-UKF Master-Slave Filter 

In this Section, the QMRPF-UKF master-slave filtering method is proposed, which uses the 
QMRPF and UKF algorithms to estimate the rotation quaternion and the gyro bias parameters, 
respectively. Section 3.1 formulates a detailed process of QMRPF master-filter. Section 3.2 describes 
the UKF slave-filter process. The flow chart of the QMRPF-UKF master-slave filter with 
multiresolution approach embedded at time/: =M is shown in Figure 1, where ® denotes the master 
filter, and ® denotes the slave filter. 

3.1. QMRPF Master-Filter 

QMRPF master-filter is provided to estimate the quaternion from pairs of vector observations, and 
computational efficiency of quaternion particle filter (QPF) is achieved by using the spatial-domain 
multiresolutional approach. This section is divided into three parts: QMRPF master-filter mathematical 
model, the complete step sum-up of QMRPF master-filter, quaternion particle filter with 
multiresolution embedded. 

A. QMRPF Master-Filter Mathematical Model 

Taking attitude quaternion as state variable, the mathematical model of QMRPF master-filter is 
built. State equation is 

q k =Q.{G) k _ l -P k _ l )q k _ 1 (11) 

Measurement equation is Equation (9). From Equation (11), it can be seen that, the gyro bias 
estimation at time k-l will be used for estimating the attitude quaternion at time k . 
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B. Complete Step Sum-up of QMRPF Master-filter 

The time is supposed to be k = M when the multiresolution approach is embedded into the 
quaternion particle filter. 

(1) Initializing the quaternion particles {q' 0 } ^ and setting quaternion weights as {w' 0 } = l/N, 
i = l,--,N. 

(2) For k<M , quaternion particle filter is used to estimate the attitude quaternion q k at time k , 
and the quaternion particles [q' k \. i and jw^j i are also obtained. 

(3) For k = M , QMRPF algorithm is adopted to determinate the rotation quaternion q k at time k , 

the number N of particles is reduced to be ./V . {q' k j ^ and jw^j ^ are derived at time M . 

(4) For k > M , the basic step is the same as step 2, but the number of quaternion particles is ./V . 



Figure 1. QMRPF-UKF master-slave filter with multiresolution approach embedded at time A; = M . 
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C. Quaternion Particle Filter with Multiresolution Embedded 
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In the following, we will show how q k is estimated when the multiresolution approach is embedded 
into the quaternion particle filter at time k=M . 

(1) Particle weights multiresolution decomposition 

A partition of N samples into N/ 2 L blocks is used where L is the level number of multiresolution 
decomposition, and then particle weights of every block are decomposed into lowpassed and 
highpassed components. 

wi, ,, 



W. 



m,lh,k-l 



w, 



-1 

m,h L ,k-\ 



m,hj ,k-\ 



W„ 



m,k-l 



m 



= l--;N/2 L 



(12) 



where T denotes an orthogonal wavelet decomposition transform matrix; m denotes the m -th block; 
\w\nk-\U=\ WQ tne particle weights of the m -th block at time k—1; w^^is the lowpass -filtered 



particle weight; \ w' mh k _A are the highpass-filtered particle weights at level j(j = !,■■■, L) , and the 



number of weights at scale j is 2 



7-1 



(2) Thresholding the highpass-filtered particle weights 

A simple thresholding is carried out on the highpass-filtered particle weights of every block: 



w. 



m,h: ,k-l 



W, 



m,hj,k-l ' 



w 



m,hj ,k- 



o, 



w. 



m,hj,k-l 



(7=l,---,L;/ = l,---,2^) 



(13) 



where T s is the wavelet threshold. Thus we have the thresholded weights: 

y±.m,lh,k-l ~ y W m,li ,/t-l ' — m,h, ,k-l ' ' "' ^m,h L ,k-l — m,h L ,k-l J 



(14) 



(3) Quaternion particles selection 

The thresholded weight corresponds to quaternion particles [q' m j htk ^) Therefore, we can 

just propagate the elements in q l mk _ x = T 1 q l miKk _ x which corresponds to the distinct elements in 
W m ,k-i = ^ T W m ,ih,k-i ■ F° r me repeated elements of W m k _ { , we can propagate one corresponding 
representative element in q l mk _ x , The number of quaternion particles is selected to be N by the 

thresholding operation, which can significantly reduce the computation in propagation while 
maintaining performance. Some detailed information about the multiresolution process can be 
referenced to [11]. 
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Now let (<Zt_i|. l and {wJUj. x denote N independent quaternion particles and their associated 
weights, ./V denotes the number of quaternion particles selected. 

(4) Time update 

The prediction values {<?L ,} are obtained from Equation (11), where f5 k x is the estimation of the 



gyro bias at time k -1 via UKF algorithm. 

(5) Measurement update (the update of weight): 



< « P hkh (h = Ps h (h -Mqi^XH-vi^-'fi (15) 

where P bk \ qt (b k q l k \ k _i) is the likelihood probability of the measurement b k associated to a given 
quaternion qL_ x , p Sbt is the probability density function of 8b k . Particle weights w{ are normalized to 



w 



(16) 



(6) Quaternion estimation 

At time k , ./V weighted quaternion samples are available. To obtain the filtered quaternion via 
minimum mean square error (MMSE) approach, the following minimization problem is solved [6]: 



II 2 

m / n Z A * ^kik-x ) - K ' subject to A k T A k = I 



'3x3 



(17) 



where A k denotes the orthogonal attitude matrix associated with the filtered quaternion, and is the 

Frobenius norm. The constrained minimization problem (17) is equivalent to the unconstrained 
maximization problem in quadratic form as follows: 



maxJr 

k 



i=i 



with: 



K = 



B k +B k ' -hJr(B k ) z 

z tr(B k ) 



max q k T Kq k 



B k (3,2)-B k (2,3) 
5,(1, 3) -5,(3,1) 
B k (2,l)-B k (l,2) 



(18) 



(19) 



Letting: 



(20) 



As in Davenport's well-known q-method, the quaternion that solves the maximization problem of 
Equation (18) is the normalized eigenvector corresponding to the largest eigenvalue of K : 
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KQk = K*A* (21) 
Then the attitude quaternion estimation q k at time k = M is obtained. 

3.2. UKF Slave-Filter 

The gyro output can't be used for attitude estimation directly for the reason that the gyro has a bias 
vector, and thus the gyro bias is considered in the estimation process. The UKF slave-filter is designed 
to estimate the gyro bias vector. Taking gyro bias as state variable, the state space model of slave filter, 
UKF, is built on the basis of gyro bias model and vector observation model. State equation is 
Equation (6), and observation equation is: 

b k = A(Q(^_! - P k _Mk-i) r k +Sb k (22) 

The noise variance matrix of state model and observation model are supposed to be Q , and R Sb , 

respectively. =cr*I 3x3 . 

It can be seen that, the state equation is linear and the observation equation is nonlinear. Moreover, 
the attitude quaternion at time k-l will be used to estimate the gyro bias at time k . The steps of gyro 
bias estimation using UKF is as follows: 

(1) Initialization 

A = *[/u p Po = £[(#> -A) (A - A) 7 ] (23) 

where f3 a and /?„ are the real value and estimation value of gyro bias at the starting time, respectively. 
P fi denotes the estimation error variance matrix at the starting time. 

(2) Calculating sigma points 

**-i(o)=A-i 

Zk-i 0" + *) = A-, - (V^T), 

where yfr , and P- denote the estimation value and estimation error variance matrix at timefc-1, 

respectively; n denotes state dimension, and for this application, n = 3 . Supposing that the estimation 
error variance matrix P- = AA T , „ IP. ) denotes the i -th column of matrix A . 

(3) Calculating the weights of sigma points 

The weights corresponding to the sigma points are given by: 

'W(O) = r/(/i + r) 

W(i) = l/[2(n + T)] ,i = \,...,n (25 ) 
W(n + i) = l/[2(n + z)] 
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(4) Time update 

The predicted mean and error variance matrix of the gyro bias at time k are given by: 

A = /t,. P k = P A _ l+ Q„, (26) 

(5) Measurement update 

V-i (?) = A(Q( ^-i " Xk-i (0)4-i K , / = 0, 1, 2n (27) 

The prediction mean of the measurement vector at time k is calculated using a weighted sum of the 
points (/) , which is given by: 



b- k = f> '(0Vi(0 



i=0 



(28) 



The measurement variance and cross correlation between /? F and £> F are computed by: 

2n 

= ZwcoCr^CO-A) (4 H (0- 4-) r (30) 

The Kalman gain is: 

K k =P Ah P h !, t (31) 
Then the estimation value and estimation error variance matrix of gyro bias at time k are described 

by: 

P k =P- k +K k {b k -b- k ) (32) 
P k =P A~ KkP hh K k (33) 

4. Simulation Results and Analysis 

In our simulation, the micro-nano satellite orbit is in sun synchronization orbit with 900 km altitude. 
Gyro and magnetometer are integrated for attitude determination. The magnetic field reference is 
modeled using a lOth-order International Geomagnetic Reference Field model. Magnetometer noise is 
modeled by a zero-mean Gaussian white-noise process with a standard deviation of 50 nT. The gyros' 
output is contaminated with measurement noise having two components: white, zero-mean Gaussian 

process with intensity ^0.31623//rad/s 5 j and gyro bias modeled as integrated Gaussian white noise 

with intensity ^3. 1623x10"* //rad/s^ . The initial attitude error angle is 0°, the initial gyro bias is 

0.2(7h) on each axis, the gyro bias variance matrix is[0.2(°/h)] 2 x/ 3x3 , the measurement noise matrix 

of magnetometer is 50 2 / 3x3 (nT) 2 , the sampling frequency is 10 Hz, simulation time is 3,000 s. The 

initial particle number of QMRPF-UKF method is 600. The threshold is 0.5x10 4 . The starting time of 
performing multiresolutional transform is the 1,000 s. QMRPF-UKF master-slave filtering is used to 
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perform the attitude determination of this micro-nano satellite. The codes are run on a computer, the 
processor speed of which is 2.5 GHz and memory is 2 Gbyte. In testing the filter, we look for two filter 
attributes: computational cost and accuracy. 

The estimation error of gyro/magnetometer integrated attitude angle estimation is shown in Figure 
2. In Figure 3, the gyro bias estimation error is given. It can be seen that the attitude angle estimation 
error is smaller than 0.1°, and the gyro bias estimation error is smaller than 0.0 17h. 



Figure 2. Estimation error of attitude angle. 
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Figure 3. Estimation error of gyro bias. 
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Here, we run QMRPF-UKF algorithm in 50 runs of Monte Carlo simulation with different 
thresholds to validate the effectiveness and complexity of QMRPF-UKF for the attitude determination 
of a micro-nano satellite. The performance is measured via RMS (root mean square) estimation error of 
Monte Carlo simulation. The smaller the RMS estimation error resulting from the filter is, the better 
the filter performance is. Figure 4 gives the selected particles numbers after performing multiresolution 
particle filter with different thresholds. RMS error and the computational cost of QMRPF-UKF with 
different thresholds are given in Table 1. It can be seen that, the larger the threshold is, the smaller the 
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number of surviving particles is. With the increase of threshold, the computational time is decreased, 
and the filtering varies in estimation precision. 

QMRPF-UKF is compared with the existing method, MLERPF, in terms of RMS error and 
computational cost, as shown in Table 1. It can be seen that, QMRPF-UKF outperforms MLERPF 
significantly in terms of computational cost and QMRPF-UKF is comparable to MLERPF in terms of 
estimation precision. For MLERPF, the state dimension of performing particle filter is six, which is 
larger than that of QMRPF-UKF. Moreover, in the calculating process of MLERPF, the particle 
number is invariant. For QMRPF-UKF master-slave filter, the filtering performance is improved from 
two aspects. On the one hand, in order to reduce the computational cost rising from the augment of the 
state dimensions, QMRPF algorithm is used as the master filter to estimate the attitude quaternion, and 
UKF is used as the slave filter to determine the gyro bias parameters. On the other hand, QMRPF 
copes with the curse of computational complexity related to the particle filtering technique by 
introducing multiresolution approach that permits a significant reduction in the number of particles. 

Three simulation scenarios are carried out to verify the convergence performance of QMRPF-UKF. 
The first one is that, the initial attitude error angle is 50°, 50°, and 50° for each axis, respectively. The 
second one is that, the initial attitude error angle is -70°, 40°, and 100° for each axis, respectively. The 
third one is that, the initial attitude error angle is -60°, 30°, and 130° for each axis, respectively. These 
are added into the initial condition of attitude estimation. For each simulation scenario, other 
parameters are same as the above mentioned ones. QMRPF-UKF convergences in each simulation. The 
attitude angle estimation error is smaller than 0.1°, and the gyro bias estimation error is smaller 
than 0.0 17h. 

Figure 4. Selected particle number after performing QMRPF-UKF with different thresholds. 
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Table 1. Comparison of QMRPF-UKF with MLERPF in terms of RMS error and 



computational cost. 





Threshold of QMRPF-UKF 


MLERPF 


1 x 10 s 


5 x 10 s 


7 x 10" 5 


1.3 x 10" 4 


1.7 x 10" 4 


1.9 x 10" 4 


5 x 10" 4 


RMS error (°) 


0.0115 


0.0137 


0.0128 


0.0108 


0.0134 


0.0124 


0.0126 


0.020 


Computational 

cost (s) 


223 


221 


219 


194 


180 


171 


160 


1840 



5. Conclusions 

We have developed a computationally efficient QMRPF-UKF master-slave filtering method, to 
solve the problem of gyro/magnetometer integrated attitude determination of a three-axis stabilized 
micro-nano satellite. Computational efficiency of the proposed filtering is achieved by combining the 
master-filter with slave-filter and by adopting the spatial-domain multiresolutional approach to reduce 
the needed number of particles while maintaining comparable estimation accuracy. The larger the 
threshold is, the fewer the number of surviving particles, and the filtering varies in estimation 
precision. 
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